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Abstract. 

It has been argued that for a finite two-dimensional classical Coulomb system of 
characteristic size R, in its conducting phase, as R — > oo the total free energy (times the 
inverse temperature (3) admits an expansion of the form: (3F = AR 2 + BR + g^mi?, 
where \ is the Euler characteristic of the manifold where the system lives. The first 
two terms represent the bulk free energy and the surface free energy respectively. The 
coefficients A and B are non-universal but the coefficient of In R is universal: it does not 
depend on the detail of the microscopic constitution of the system (particle densities, 
temperature, etc.). By doing the usual Legendre transform this universal finite-size 
correction is also present in the grand potential. The explicit form of the expansion 
has been checked for some exactly solvable models for a special value of the coulombic 
coupling. In this paper we present a method to obtain these finite-size corrections 
at the Debye-Hiickel regime. It is based on the sine-Gordon field theory to find an 
expression for the grand canonical partition function in terms of the spectrum of the 
Laplace operator. As examples we find explicitly the grand potential expansion for a 
Coulomb system confined in a disk and in an annulus with ideal conductor walls. 



PACS numbers: 05.20.Jj, 51.30.+i 
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1. Introduction 

There are several reasons to study models over a finite-size region. For instance, with 
the recent advance of computers, much information on statistical models has been 
derived from computer simulations, which are necessarily limited to systems of finite- 
size. Also experimental systems are finite (although very large). The finite scaling 
hypothesis allows the study of some response functions for such finite-size models. For 
a ci-dimensional system, this finite scaling hypothesis states that if a given response 
function (for example the susceptibility in a magnetic system) diverges in the bulk like 
£ d ~ 2x , where £ is the correlation length and x is called the scaling dimension of the 
corresponding studied quantity (the magnetization in the above example), then in a 
finite system of characteristic size R, the response function should obey the scaling law 
R d ~ 2x Q(R/l;), where $ is some universal scaling function PflEj- At the critical point 
where the correlation length diverges, the response function is then proportional to 
R d ~ 2x . However the scaling of the free energy at criticality is less well understood at 
least for arbitrary dimension. In two dimensions, using methods from conformal field 
theory, it has been shown that for a finite system with smooth boundary, of characteristic 
size R as R — ► oo, at criticality, the total free energy F has a large- R expansion of the 
form [3 IS] 

/3F = AR 2 + BR- ^lnR + 0(l) . (1) 
6 

with = (keT)- 1 the inverse reduced temperature. The first two terms represent 
respectively the bulk free energy and the "surface" (perimeter in two dimensions) 
contribution to the free energy. In general, the coefficients A and B are non universal 
but the dimensionless coefficient of In R is universal depending only on c, the conformal 
anomaly number, and on Xi the Euler characteristic of the manifold (x = 2 — 2h — b, 
where h is the number of handles and b is the number of boundaries). Surprisingly 
enough, for classical Coulomb systems in their conducting phase — not at criticality — 
this expansion for the free energy seems to holds with c = 1 and a change of sign in 
the last term, and it has been explicitly checked for Coulomb systems lying on some 
simple geometries for some exactly solvable models for the fixed temperature defined 
by (3q 2 = 2 where [3~ l = ksT and ±q are the charges of the particles |U El El E| and 
also it has been verified numerically for the one-component plasma confined in a disk jH] 
for other values of the coulombic coupling. The existence of the universal finite-size 
correction has also been proved for the two-component plasma confined in a disk with 
hard walls in the whole regime where the system of point particles is stable (/3q 2 < 2) 9J. 
For the one-component plasma JU], for the symmetric two-component plasma [TJ and 
for the asymmetric two-component plasma ^2] confined in a sphere the existence of the 
finite-size correction has been proved for any temperature (provided that the system is 
stable) by application of the stereographic projection and some non-trivial sum rules 
concerning the density-density correlation function in the plane geometry [T3| I14j. 
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Although the particle and charge correlations in Coulomb systems in their 
conducting phase are short-ranged because of the screening, the electric potential 
correlations are long ranged [13 EE!- It has been argued [H E] that in this sense 
the system is comparable to a critical system and therefore the expansion of the free 
energy (JTJ) should hold. 

In this paper we present a method to obtain the grand potential finite-size expansion 
for a confined Coulomb system in the Debye-Hiickel regime and we check the existence 
of the logarithmic universal finite-size correction. The Debye-Hiickel regime is defined 
by the requirement that the average coulombic energy is much smaller than the thermal 
energy ^3 EE] ■ By the usual Legendre transform between the free energy and the grand 
potential it can be inferred that the free energy will have the same logarithmic finite- 
size corrections as the grand potential. Our method is based on the sine-Gordon field 
theory ^H] to calculate the grand canonical partition function. 



2. Sine-Gordon theory in the Debye-Hiickel regime 

There is a well-known analogy between statistical mechanics and quantum field theory: 
often the partition function of a d- dimensional statistical model is formally analogous to 
the generating functional of a quantum field in d space-time dimensions in the Euclidean 
formalism [20] • The simplest example of a quantum field theory which has a relevance in 
statistical mechanics is the free Boson or Gaussian model. In this section we show that 
the grand canonical partition function of a Coulomb system in the Debye-Hiickel regime 
may be represented as the generating functional of a massive free Boson theory and 
therefore it can be obtained exactly from a Gaussian integration as an infinite product 
of functions of the eigenvalues of the Laplace operator calculated on the manifold where 
the system lives. 

Let a classical (i. e. non-quantum) Coulomb system be composed of a = 1, . . . ,r 
species of particles each of which have N a charges q a confined in a ^-dimensional 
Riemannian manifold of volume V. Suppose that the system is confined by grounded 
ideal conductor walls, thus imposing vanishing Dirichlet boundary conditions to the 
electric potential. We shall describe the system in the grand canonical ensemble with 
fugacities ( a = e^ Ma /A d for the species a, where fi a is the chemical potential and A is 
the de Broglie thermal wavelength which appears when the Gaussian integration over 
the kinetic part of the hamiltonian is carried out. For a finite but macroscopically large 
system, the interior of the system will be at an almost constant electric potential ipo- 
The value of ip is controlled by the choice of the fugacities. We will suppose in the 
following that the fugacities satisfy the relation 

= (2) 



which is often referred in the literature [2H] as the pseudo-neutrality condition. In 
the | Appendix B we consider the general case when the fugacities do not satisfy the 



condition (J2J). 
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Let us introduce the Coulomb potential for a non-confined system for unit charges 

-^—7 , if d = 3 



OH |r ~, ' (3) 

-In 1 - 1 , if d = 2. 

In two dimensions, a solution of Poisson equation that vanishes at large distances does 
not exist, therefore it is necessary to introduce an arbitrary length L that fixes the zero of 
the electric potential. However as we will see later, in the formulation of Debye-Hiickel 
theory proposed here it will be necessary to suppose that L — > oo, thus receding the 
zero of the electric potential to infinity. The necessity to choose L — > oo also appears 
in the formulation of Debye-Hiickel theory from Ornstein-Zernicke equation and the 
approximation of the direct correlation function by the Coulomb potential [2*Tj . 

The Coulomb potential in d dimensions for the system confined with Dirichlet 
boundary conditions will be referred as Vd- This potential satisfy the Poisson equation 

Av d {r,r') = -s d 5{r-r') (4) 

with S2 = 2tt and S3 = Air, and the Dirichlet boundary condition. If one considers 
v d (r, r') as the kernel of an integral operator which we will also call v d , we have 
Vd = — s^A -1 . Let \I / ri (r) be the normalized eigenfunctions of the Laplacian with 
Dirichlet boundary conditions, that is A\I/ n (r) = A n \P n (r) where A„ < are the 
corresponding eigenvalues. These functions are also eigenfunctions of v d with the 
corresponding eigenvalues —Sd/X n > 0. Consider two particles located at v a ^ and rpj. 
A standard operator spectral decomposition gives for the interparticle potential 

Wrf(r a ,i, Tpj) = - V" *n(ra,i)*n(r^j) • (5) 
— A n 

n 

The bar over \I/ indicates complex conjugation. Additionally to the interparticle energy 
we consider the energy of each particle located at r a j in presence of the field produced 
by itself i^(r a) j, r a) j) = Vs-E( r a,i)- This term is (twice) the self-energy of a unit charge. 
Proceeding as in (JSJ it may be given by 

vs- E (r a ,) = -Y,^K( r ^\ 2 ( 6 ) 

k X k 

here the A° refer to eigenvalues calculated for the system without boundaries. Of course, 
because of the form of the Coulomb potential the self-energy is in fact infinite. This 
divergence may be avoided by the introduction of a short distance potential [TH] to 
cutoff the singularity of the Coulomb potential at the origin. To simplify the notation, 
we will not write down explicitly this short distance potential in what follows. It should 
be noted however that the introduction of this short distance potential is mandatory 
in classical statistical mechanics of Coulomb systems in order to have a well defined 
partition function in three dimensions (in two dimensions at low coulombic couplings 
a system of point particles is stable). On the other hand it turns out that the Debye- 
Hiickel approximation is well defined for a system of point particles: as we shall see later 
the short distance potential does not appear in the final results. Let us remark that 
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for systems governed by quantum mechanics if all particles of a same sign are fermions 
then the system is stable [22] • This is the case in nature where quantum mechanics are 
responsible for the creation of stable bound states. Therefore our classical analysis will 
apply only to fully ionized systems. 

The potential energy of the system is given by 

1 x ' I 

H = o Z~Z 9aQpVd(r a ,i, rpj) + - V V ql [v d (r aji , r a>i ) - Us_jg(r tt)i )] . (7) 

Z h3 Z 

a,p a i 

The prime in the first summation mean that the case when a = (3 and i — j must be 
omitted. The first term is the usual interparticle energy between pairs. The second 
term is the Coulomb energy of a particle and the polarization surface charge density 
that the particle has induced in the boundaries of the system. When the method of 
images is applicable to compute the Coulomb potential Vd, this energy can also be seen 
as the energy between each particle and its image. 
Using the microscopic charge density defined as 

r N a 
a=l i=l 

we can write the potential part of the Hamiltonian of the system as 

Iff 1 r Na 

H =o dr dr' p(r)v d (v,v')p(r') - - ^ g^^fr^) . (9) 

•* V J V a=l i=l 

Notice that with this notation, in the first term, the terms q^Vd(r ati ,r a ^)/2 have been 
included. Often in the sine-Gordon transformation 123] the second term is omitted 
in the Hamiltonian in equation Q, which implies that the self-energy (infinite for point 
particles) is included in the Hamiltonian when it should not be. This problem can be 
cured with a proper renormalization of the fugacity [21] , however this method is not the 
more convenient to use for the Debye-Hiickel approximation, therefore we will proceed 
to subtract the self-energy from the start as shown in equation J5J). 
Now, using the well-known Gaussian integral 

5 B-A -B = J r 1Q \ 

JdXe-^ x - A - x 1 ^ 

we can represent the Boltzmann factor asf 

N, 



e ^ H = < exp 



/R 
zp(r)0(r) dr + - ^ £ v s-e(* 
a=l i=l 



where we have defined the average of any operator 6 as (o) = 4 j T>(f)de 2s <i /^( r ) A ^( r ) dr ; 

with Z = J V(j) e 2s d /^( r ) A< ^( r ) dr _ Q n the other side, using (jHJ) and (fTTj) after some 

| Rigorously speaking this Gaussian transformation can not be formulated with the Coulomb potential 
Vd(r, r') because this potential diverges when r = r'. This problem can be solved as in Ref. using 
instead a potential like (1 — e~ Kr l e )jr which is regularized at short distances and taking the limit e — > 
at the end of the calculations. Again, for simplicity, we will omit explicitly this detail in what follows. 
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calculations the grand partition function is given by jTH| 



1) 



Nx\...N r 



N r 



JVi=0 



cxp 



iV r =0 



a=l j=l 



E 



-iq a 



dr 



(12) 



with V the volume of the manifold containing the system. We see from equation (|12j) 
that the partition function for a gas of particles with two-body interactions may be 
obtained as the average of the partition function of an ideal gas in an external fluctuating 
potential i<p(r). In general the integration involving the calculation of ()12|) cannot 
be performed analytically — with the notable exceptions of the two-dimensional two- 
component plasma (symmetric 1 : 1 and asymmetric +2 : —1) which has been exactly 
solved in the bulk and in some semi- infinite geometries [2U 123 I2M ■ 

The coulombic coupling T is defined as the average Coulomb energy divided by 
the thermal energy. We can actually define a coupling for each species of particles as 
follows. In two dimensions these are defined as T2. a = f3q&- On the other hand, in three 
dimensions T 3 a = Pq 2 (a 3 . In the Debye-Hiickel regime we have T djCt <C 1 and we can 
expand 



exp 







q 2 



E {r) = l~t(3q a (j){r) +—vs-E{r) +o {T da ) .(13) 

In two dimensions the field 0(r) has dimensions of charge q a therefore the above 
approximation is an expansion to the order (r 2 , a ) 2 in the coulombic couplings r 2)Q . In 
three dimensions the field 0(r) has dimensions of charge / distance. One can do a change 



of variable in the functional integral to have a dimensionless field 



c 



-1 



q, then 



it is clear that the approximation (fT3*|) is again an expansion to the second order in the 
coulombic couplings T 3a . Notice that the self-energy term /3q 2 vs-E( T ) is already of 
order (T^a) 2 - This can be seen by noticing that the covariance of the Gaussian measure 



W(j)e 



is (0(r)0(r')> = (3- 



_1 v d (r, r') 



Then it is clear that the self-energy term 



/3g 2 t>5_£(r)/2 is of the same order as (/3g Q 0(r)) 2 . Therefore we do not include any 
terms of order (/3g 2 f U5_ E (r)/2) 2 or p 2 ql i (p(r)vs-E( r ) which are of higher order in the 
expansion fHTj) . 

Using equation (|T3*|) and the pseudo-neutrality condition (j2J) we have 



E(/3,Ci-C*,^) = (exp 



(a ((3q a (j)(r)) 2 
9 dT 



L^i a — 1 Z-j i 



2,\ 



TT 



(14) 



where the spectral decomposition (jH)) of the self-energy and the normalization condition 



dr — 1 have been used to write the contribution the self-energy terms as a 



sum over the eigenvalues A° of the Laplacian without boundaries. Now performing the 
Gaussian integration the averaged quantity equals to 



1 



V(j)e z 



i/0(r)(M-E a C Q (/3fc) 2 )*(r) 



dr 



det 



1 - 



Ea S dCaPql 



-1/2 



(15) 
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Notice that the averaged quantity we had just computed is the generating functional 
of a free Boson theory with mass proportional to the inverse Debye length defined by 
k = \/J2 a s d(,af3<la- Using the invariance of the determinant we obtain for the grand 
canonical partition function 

s(/3,c 1 ,...,c,^)= [n(i-£)ri e *) 1 eZaV(a - 

This is simply a product of factors that are a function of the eigenvalues A, of A. The A^s 
depend on the shape of the domain in which the Coulomb system lies. As we see, they 
constitute a natural way to introduce the information on the domain to calculate the 
corresponding finite-size expansion of the grand potential Q. It is interesting to point 
out here that in the case of a non-confined system A n = A° and the infinite products 
in (|16j) become a regularized Weierstrass product Y[ n fl — f?rj e K2 / A ™ . The Yl e K2 / A ™ 

term cancels out the ultraviolet divergence coming from of the Y[ fl — y~j term. This 
product converges for systems in three dimensions (301 • However, as it will be seen in the 
next section, in two dimensions some infrared divergences appear and the product must 
be regularized by introducing a lower cutoff. The sine-Gordon transformation has been 
known for some time [T%1 [T§\ |2"3|. For three dimensional non-confined systems the sine- 
Gordon transformation have been used to go beyond the Debye-Hiickel approximation 
and to perform low fugacity [31] , high temperature [33] or loopwise j3H] expansions. The 
main point of this section was to show that the proper subtraction of the self-energy 
terms (which have to be added initially to perform the sine-Gordon transformation) leads 
to a well-defined expression for the grand potential in the Debye-Hiickel approximation 
which could be eventually evaluated for confined systems. In the |Appendix A we show 
how this formulation of Debye-Hiickel theory is related to the usual one. In the following 
section we apply this method to the calculation of the grand potential Q = — A^TlnH 
of a Coulomb system confined in some simple geometries. 



3. Finite-size corrections to the grand potential for a confined Coulomb 
system 



3.1. N on- confined systems: the bulk 

Before applying our method to confined systems let us illustrate some points of the 
calculation of the grand potential from equation (jlfij) for a bulk system. For a d- 
dimensional non-confined system the wave functions corresponding to the A° are given by 



^ and A*°(r) = A k 
the grand potential is given by 

2 J (2nY 



-kX(r) 



Then replacing into equation (Ti7j|) 



In 1 + 



K 



K 



21 



(17) 
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In three dimensions d = 3 the above integral is convergent giving the result for the bulk 
grand potential per unit volume 



8fl K 3 > 

1t = -isf-5> < 18 > 

a 

The density n a of the particles can be obtained by using the usual thermodynamic 
relation 

n Q = -<«^P- = Ca + (19) 

which can be replaced back into equation (|18|) to give the well-known equation of state 
from Debye-Hiickel theory [23 EH] 

Notice that in the last equation the Debye length that we have defined by k^ 1 = 
E Q s d(af3qa] ^ 2 has been replaced by the usual Debye length in terms of the density 
k dh = [Ea s d n a/3qa] ^ 2 ■ This is correct at the order of approximation we are working 



3/2^ 
Z,a> 



on since k = «dh[1 + 0(T 

Let us point out that the proper substraction of the self-energy terms makes the 
integral (fTTj) convergent and avoids the need to use some other arbitrary regularization 
scheme, like for instance dimensional regularization used in Ref. [2H] which by the way 
yield the incorrect result J2 a Ca — k 3 ' j (247r) for the pressure when it is expressed in 
terms of the fugacities. Our regularization scheme is actually equivalent to the normal 
ordering for the product : 0(r) 2 : used in Ref. [26J. 

As it was pointed out in the preceding section the infinite product (|TH|) for a non- 
confined system is a regularized Weierstrass product. The order of the sequence of the 
Laplacian eigenvalues is // = d/2 [HO], therefore for d = 3, ji = 3/2 > 1 with integer part 
equal to one, and the terms exp(/t 2 /A°) in the product (fTtjj) are enough to regularize the 
infinite product. 

The situation in two dimensions d — 2 is more delicate since jj, — 1 is a limiting 
case. If we blindly try to compute (JT7j) we will notice that the integral is not well-defined 
for k — > 0. Trying to cure an ultraviolet divergence, we introduced an infrared one. The 
problem can be traced back to the spectral decomposition (J3J of the Coulomb potential. 
Evaluating explicitly the interparticle energy using the expression (j^J) gives 

i p2tt poo j ifc|r-r'|cos6» poo jj, 

v 2 (t,t>) = ^ ^ ^ dkdB = J y/o^ir-r'l) (21) 

with Jo the Bessel function of order 0. This integral diverges at k = 0. To avoid this, 
we introduce a cutoff k min at k — > 

v 2 (r,r')= ^ ^J (k\v-r'\) (22) 

^min 

= -C + ln . _\ +o(l) (23) 
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where A; min — > and C is the Euler constant. Since we know that t>2(r, r') = — In ^ * ! 



we can find the expression for fc n 
— In ^ r ~ r ■ ; then 

2e~ c 



by comparison: i^(|r — r' 



In 



lr-r'|fc m 
2e-C 



L 



(24) 



We note that for the above calculation being consistent it is necessary to choose L — > oo, 
as it has been discussed earlier in the preceding section. The necessity of this choice for 
the arbitrary constant L is also discussed the appendix B of Ref. |21| . 

Returning to the calculation of the grand potential in two dimensions, we impose 
the infrared cutoff k m i n = to the integral (jT7jl to obtain the result for the grand 
potential 



/3n 



K 

Air 



, kL _ 1 

-In C + - 

2 2 



In the above expression all terms that vanish when L 
density-fugacity relation is now 



(25) 

oo have been omitted. The 



n. 



-c 



d(pn/v) 



C«-C 



Hi 



(26) 



d( a w 2 

For a two-component plasma it can be checked that this result is reproduced from 
the small- (3q 2 expansion of the exact relation between the density and the fugacity |24j . 
Notice again that k can be replaced by kdh at the order of approximation we are working 
on, since k = kdh[1 + 0(T2, a In r 2 , a )]. Reporting equation back into equation 
one obtains the equation of state, which turns out to be exact at the level of the Debye- 
Hiickel approximation, 



pp = Ua ( 



1 - 



(27) 



Doing the usual Legendre transform F = Q + ^Z a ^ a N a , one can recover the known 
expression for the excess free energy in the Debye-Hiickel approximation 



PF e: 



K 



DH 



1 , K DH-^ ^ 

In C 

2 2 



(28) 



V 4tt 

To conclude with the results for a two-dimensional system let us clarify a point 
regarding the limit L — > oo. Actually in equation ()23|) and below we require that L 
be large compared to the average distance between particles which is of order n~ 1 / 2 
with n the density. In the Debye-Hiickel approximation the density n is of the same 
order as the fugacity £. Therefore we require that LC 1 ! 2 ^> 1. In the results for the 
grand-potential (|23j) . the densities (|2*B^1 and the free energy pSjl appears the quantity 
kL which is proportional to (f3q 2 ) 1 ^ 2 L( 1 ^ 2 = r 1 / 2 (L^ 1 ^ 2 ). Notice that in the above 
expression LC 1 ! 2 ^> 1 but the coulombic coupling r< 1. Therefore we require that in 
two dimensions the Debye-Hiickel limit should be taken with T — > 0, LC 1 ! 2 — > oo but 
r 1 / 2 (L(^ 1 / 2 ) should remain of order 1. 
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3.2. The disk 

We now consider a two-dimensional Coulomb fluid confined in a disk of radius R. To 
apply the method outlined in section we first need to compute the eigenvalues of the 
Laplace operator for this geometry. Let ^(r,(p) = R(r)<&(ip), we look for the solution 
of the equation A^f(r,(p) = X^f(r,tp). Using the explicit form of the Laplace operator 
in polar coordinates we find ^f(r,(p) = R(r)<f>((p) oc e ±llv Ii{\f\r) where h(x) is the 
modified Bessel function of order I. Using the boundary conditions *&(R, (p) = 0; 
*(r,0) = *(r, 2ir), we find I G Z and h (\/Afci?) = 0, that is a/A^-R = u i>n is the ra- 
th zero of ij.§ Then replacing these eigenvalues into equation ()16j) gives for the grand 
potential the expression 

Z = -00 n=l \ '' n / V 7 ■'Kmin Q 

The second term, written as an integral over k, comes from the terms involving the 
Laplacian eigenvalues for a non-confined system: e K2//A °, with the same infrared cutoff 
fc min discussed previously and given by equation (J23j) . Both the sum and the integral 
diverge separately for large values of I and k but when we put together both terms 
the ultraviolet divergences should cancel. It is however more convenient to compute 
separately each term. Therefore we will impose an upper cutoff N for I and k max for k. 
Both cutoffs are of course proportional, the exact relation between N and fc max can be 
obtained at the end of the calculations by imposing that for the bulk grand potential 
we should recover the result fl25(l from last section. 

Using the infinite product representation of the modified Bessel function 
]l„=i (l - ^f) = *! [37 and the property I t {x) = U{x) we have 

N , v TV N 

pn= VlnZ! + ln — y"Z + y)lnJ,(Ki2)--lii[Jo(«i2)] (30) 

1=0 v 7 1=0 1=0 

Using Stirling approximation: In TV! = TV In TV — TV + |ln(27rTV) + • • Euler-McLaurin 

summation formula: ^Zo /(0 = So /(0<« + I [/(°) + /(# )] + ^ [/'(#) - /'(0)] + ■ ■ ■ 
and the uniform Debye expansion |HZj for In Ii(z), valid for large z, 

^Ii(z) = ~\ ln(27r) - \ In (z 2 + Z 2 ) + r,(l, z) + + o (32) 

,(/,z) = (z 2 + / 2 ) 1/2 -/sinh- 1 Q; «= f 1/a (33) 
after some calculations we finally obtain the large- i? expansion 

§ Notice that since the zeros of // are imaginary then y/Xk is imaginary, this is expected since the 
Laplacian eigenvalues A^ are negative. 
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In the above expression all terms that vanish when N — > oo and /c max — > oo have been 
omitted. The bulk term (proportional to ttR 2 ) of the above equation (}3"4"|) should be 
the same as in equation (J2H)) therefore the cutoffs N and k mSuX should be related by 
*W = f e 1 ' 2 - Finally 

pQ = [3u b 7cR 2 + + - MkR) + 0(R°) (35) 

6 

with the bulk grand potential per unit volume u b (equal to minus the bulk pressure p b ) 
given by 



K 2 



(3uj b = -(3p b 

47T 



- In C + - 

2 2 



and the surface tension 7 is given by 

Pi = ~ ■ (37) 

The two-dimensional two-component plasma near a plane ideal conductor electrode has 
been solved exactly |2Zj. For a two-component plasma our result ()37)1 for the surface 
tension agrees with the lower order expansion in f3q 2 of the exact result of Ref . • In 
equation (|35jl we notice the existence of the universal logarithmic finite-size correction 
(1/6) In R with \ = 1 for the disk. 

3.3. The annulus 

We now consider a Coulomb fluid confined in an annulus of inner radius R\ and outer 
radius R 2 . As before we need to calculate the eigenvalues of the Laplace operator for 
this geometry. The eigenfunction of the Laplacian with eigenvalue A, in this geometry, 
is ty(r,(p) = AI[(\/~\r) + BKi(\f\r) e llLp . Imposing the Dirichlet boundary conditions 
yields the linear system of equations *ff(Ri,(p) = and ^(R 2 ,(p) = 0. To have a non- 
vanishing solution for the eigenproblem we require that the determinant of this system 
vanishes. This gives the equation that defines the eigenvalues for this problem, 

I l (V\Ri)K l (y/XR 2 ) - K l (y/\R 1 )I l (V\R 2 ) = (38) 

this means that = zf n where z^ n is the n-th root of equation 

I^zR^K^zR,) - K l {zR 1 )I l (zR 2 ) = . (39) 

Notice that the roots of this equation are the same for I and — I, therefore we will 
concentrate on the case I > 0. To compute the grand partition function from 
equation (fTtjj) we need to evaluate the infinite product Yli Tin ~ «f~)" "^ or a S^ ven 
I, the product over the index n of the roots of equation (|3*9*j) can be performed using a 
generalization of the infinite product representation of the Bessel functions used in the 
last section |31[7|EH1- For I > 0, let us introduce the entire function 

M*) = , ^ , x 1 [IiizRiWzRz) - Ki{zRi)Ii(zR 2 )} . (40) 
Sx 1 _ [ Ei 



R2) \R 
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By construction the zeros of the function /; are z^ n and it has the following properties: 
fi(0) = 1, //(0) = and fi(z) = f\(—z). Therefore fi admits a Weierstrass infinite 
product representation of the form [39 

/'(*) = II f 1 - ir ) ■ (41) 



Z Ln 



Then the infinite product we wish to evaluate is simply Yl n (^ ~ K<2 / Z fn) = fi( K )- F° r 
I = the function f should read 

= i (v/n S [h{zRi)K {zR 2 ) - K (zR 1 )I (zR 2 )} . (42) 

The grand potential is then given by 

^ = ^ln/ ; (^) + 2 ln /o(«:)-^f / y -E^ Q . (43) 

Z = l ^ fc min a 

As in the case of the disk we regularize the summation on / by introducing an upper 
cutoff iV and the integral with an ultraviolet cutoff fc max . These cutoffs are proportional 
in order to cancel the divergences. However their exact relationship is a priori different 
from the one in the disk case and can be found at the end of the calculations by requiring 
that we recover the same bulk value of the grand potential as in the previous examples. 



On the other hand the infrared cutoff k min = ^j— is the same as before. 

We now proceed to find the finite-size expansion of the grand potential. We consider 
a very large annulus with R 1 — >• oo, R 2 — > oo, R 2 — R\ — > oo and x — R\j R 2 < 1 finite 
and fixed. The calculations are similar to those of the disk, we now use the uniform 
Debye expansion of In Ki(z) valid for large arguments [37j 



In KAz) = In 



V2 



hn{l 2 + z 2 )- V {l,z)+ln 



3m — 5m 3 



24/ 



-oItt^) (ID 



P + z 



with r)(l,z) and u defined in equation (J3~3~j) . Notice that in the functions fi(n), the 
contribution of Ki(k,R 2 )Ii(kRi) is exponentially smaller than the one from the term 
Ii{kR 2 )Ki{hRi) . Using again the Euler-McLaurin summation formula to transform the 
sum over I into an integral, after some calculations we find in the limit iV — > oo, 

(3Q = \ (r\ In -™- - R\ In -^-) (45) 

.Mg^).^^^. (46) 

All terms that vanish when N — > oo have been omitted. To recover the proper bulk 
value of the grand potential and ensure extensivity, the first term in equation (J43|) should 
vanish. This imposes the relationship between the ultraviolet cutoffs iV and k max , 

2Ne 1/2 * 2 

= R 2X ^ . (47) 

""max 

This relation is similar to the one found in the disk replacing R by R 2 x^- X . 
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Returning to the grand potential we conclude that its the large annulus expansion 

is 

tt = 7i(R 2 2 - R\)uj b + 2tt( j Ri + #2)7 + Oil) (48) 

with the Mb and 7 given by equations (jHEj) and (jHZj) respectively. In the 0(1) neglected 
terms there are terms of the form \n{Ri/ R 2 ) and more generally functions of x = R1/R2 
which are indeed of order 1. There are not logarithmic finite-size corrections, such as 
ln(K\/ R1R2) , according to the fact that \ = f° r an annulus. 



4. Summary and conclusion 

The method presented here gives a practical prescription for the calculation of finite-size 
corrections of the grand potential of a Coulomb system in the Debye-Hiickel regime, 
that can be easily applied to more complicated geometries in two and three dimensions. 
The proper substraction of the self-energies avoids the divergence of the infinite products 
involved in the calculations. In the disk and annulus geometry that we used to illustrate 
our method, we recovered the bulk pressure and the surface tension of the system in the 
Debye-Hiickel regime. For the disk we obtained a universal finite-size correction | In R, 
with the expected value x = 1> f° r the Euler characteristic of the disk. For the annulus 
since \ = no finite-size correction is expected and we confirmed this result by direct 
calculation of the finite-size expansion. In the case of a system in a domain of arbitrary 
shape, the logarithmic universal correction to the grand potential may be obtained from 
the asymptotic properties of the spectrum of the Laplace operator and its relation with 
the geometry of the manifold for which this spectrum is calculated jJOJ • Work in 
this direction is in progress. 
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Appendix A. Relationship with the usual formulation of Debye-Hiickel 
theory 

The usual formulation of Debye-Hiickel theory [3T], for a confined Coulomb system with 
Dirichlet boundary conditions for the electric potential, starts by computing the electric 
potential $ a ( r , r ') created at r' by a particle of charge q a located at r and its polarization 
cloud. We have <£> Q (r, r') = q a K{r, r') with the Debye-Hiickel kernel K that satisfies 

(A-4 H )ir(r,r') = -^(r-r / ) (A.l) 
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with 



f3q^n a Sd- Formally K can be written as 



K{r,v') 



A 



"DH 



(A.2) 



where the Laplacian is considered to satisfy Dirichlet boundary conditions. Then, one 
computes the internal potential energy U of the system as 



which can formally be written as 



Va r, r 



0) 



K 



DH 



2/9 

^DH 

2/3 



Tr 



E 



q a n a s d 



A 



K 2 



c/r 



A «|, H 



1 

A° 



1 " 

A° 
1 

~ AH 



(A.3) 

(A.4) 

(A.5) 
(A.6) 



A n n, DH 

The notation A denotes the Laplacian operator with free boundary conditions. 

On the other hand the internal excess energy U can be computed from the 
thermodynamic relation U = — (91nS/9/3)^ v . Using the sine-Gordon formulation, 
we can obtain an independent expression for the internal excess energy and compare it 
to equation (|A.6|) . Using equation (fTHjl gives 



9 In H 

dp 



2^5 



E 



In 1 - 



K 



.2 



Yp 



E 



A, 



An 
1 

" A2 



+ 



(A.7) 



At the Debye-Hiickel level of approximation kdh (expressed in terms of the densities) 
can be replaced by k (expressed in term of the fugacities) with corrections of higher 
order. Therefore with equation (jTUj) for the grand potential and equation (|A.7|) . one 
recovers the expression (jA.6|) for the internal excess energy obtained from the usual 
formulation of Debye-Hiickel theory. 



Appendix B. On the pseudo-neutrality condition and the potential 
difference between the system and the walls 

Coulomb systems have the interesting property that any excess charge in the system 
is expelled to the boundaries [221 • Therefore any infinite system is neutral. When 
the system is described in the grand canonical ensemble with fugacities the 
electroneutrality has the consequence that the fugacities are not independent. Several 
choices of the fugacities can describe the same system. More precisely, the grand 
potential does not depend on the combination ^2 a q a Ca EH1 1221 H2j- Therefore one 
can impose the so-called pseudo-neutrality condition 

Z>«£ = °- (B.1) 
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For a confined Coulomb system the situation is more involved. Suppose that the 
confined system, described in the grand canonical ensemble with fugacities ( a , is in 
equilibrium with an infinite neutral reservoir at zero electric potential with fugacities 
C* that satisfy the pseudo- neutrality condition (|B.1|) . Let us consider that the confined 
system is large and that far from the boundaries the average electric potential of the 
system is a constant tpQ. Writing down the equilibrium condition that the electrochemical 
potentials of the system and the reservoir should be equal yields C* = C <a e~ l3qo '^ . 
Therefore the confined system can be described with the fugacities ( a which a priori 
do not satisfy the pseudo-neutrality condition or with the fugacities which satisfy 
the pseudo-neutrality condition plus the parameter ipo which is the potential difference 
between the system and the reservoir. In this article we have supposed so far that the 
fugacities ( a satisfy the pseudo-neutrality conditions. In the following we will consider 
the general case when the fugacities do not satisfy the pseudo-neutrality condition and 
we will explore how they are related to the potential i/jq ^ in an approximate mean 
field picture. If ipo ^ this potential difference will create a surface charge density near 
the boundaries and we would expect that this effect will add to the grand potential 
and the free energy a surface term. We will show that this contribution turn out to be 
— (l/2)Qip where Q is the excess charge of the system which is spread over the surface 
of the boundaries 

Actually we can justify this argument within our formalism of the sine-Gordon 
transformation by adapting some arguments put forward in Ref. J2H| for the case of free 
boundary conditions. In the case of two dimensional systems, we will also show that 
if the potential difference ipo is not too high, the contributed surface term is of higher 
order in the coulombic coupling constant than the surface tension already computed in 
section 13.21 and given by equation (J3*7)l and therefore it can be neglected in the Debye- 
Hiickel approximation. 

Let us rewrite equation (|12j) for the grand canonical partition function as 

S(/U a ,V) = ^ / V<t>exp[-S] (B.2) 
with the action S given by 



S 



2s d 



(B.3) 



To lighten the notation, in this appendix we will often omit the variable r in the 
integrals: J 4> = J 0(r) dr. For simplicity we have omitted the self-energy term which is 
irrelevant in the present discussion (one could eventually consider that the fugacities ( a 
are renormalized by a multiplication by e /3q °' vs - E ^ 2 ). 

If the fugacities do not satisfy the pseudo-neutrality condition, the stationary point 
of the action S is not = as before. Let ip( r ) be i times the solution of 5S/5(p = 0. 
This field satisfies 

AV(r) + s d C*q a e- M{r) = . (B.4) 
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with Dirichlet boundary conditions: ^(r) vanishes on the boundary. This is Poisson- 
Boltzmann equation and the field ip(r) is the average electrostatic potential in the mean 
field approximation j33j- Notice that if the fugacities ( a satisfy the pseudo- neutrality 
condition then ^(r) = is a solution of Poisson-Boltzmann equation (|B.4|) . In this case, 
and in the mean field approximation, the potential difference between the boundaries 
and the system is zero. If the fugacities ( a do not satisfy the pseudo-neutrality condition 
then ip(r) is not zero and contrary to what has been done before the expansion to the 
quadratic order of the action S should now be done around = —iip instead of = 0. 
To accomplish this let us do the change of variable in the functional integral 0' = <f) + iip- 
We have T>cf) = T><p' and the action is now given by 

-S= f — [0'A0' - 2z0'A^ - ipAip) + C a e- pq ^ +l ^ . (B.5) 
J 2s d 

The new field 0' fluctuates around and now we can expand the exponential to the 
second order in the coulombic coupling. The linear terms in 0' in the action S are 
canceled by applying the stationary condition (Poisson-Boltzmann) equation (jB.4|) and 
we find S = Si + S 2 + o(T d J^) with 

Si = lj 0'(r) + (Pq a ) 2 (*e-^A 0'(r) dr (B.6) 

dr. (B.7) 



S 



A^ (r)A ^(r)-^C Q e-^ r ) 



The term Si is of order T d J^ in the coupling constants. To verify this, note that in two 
dimensions the field 0' can be written as qf(nr) with / some function of order one and 
q is the magnitude of the elementary charges in the system, for example q = max \ q a \. 
Rescaling the distances in the integral by the inverse Debye length k shows that Si 
is of order T^a- in three dimensions 0' = qnffar) and doing the same scaling in the 
integral as above shows that Si is now of order T 3 a . To know the order of magnitude 
of 5*2 we need further assumptions. To proceed, we shall need in principle the solution 
ip(r) of Poisson-Boltzmann equation (jB.4|) . However the solution of Poisson-Boltzmann 
equation is not known explicitly except for a few very simple geometries E3 | l47) | 14 7j . 
Nevertheless the qualitative behavior of the mean field ip(r) is very simple. It vanishes 
on the boundary and a few screening lengths away from the boundary it is almost equal 
to a constant value ip$. This constant average value of the potential ipo is given by 
Poisson-Boltzmann equation (|B.4|) for a constant field: 



X)feC«e~^ = 0. (B.8) 



Let us define the renormalized fugacities = ( a e~ l3Qa ^°. By equation (jB.8|) these new 
fugacities satisfy the pseudo- neutrality condition (JB.lj) . The physical interpretation of 
these new fugacities is of course the one exposed at the beginning of this appendix: 
they are the fugacities of the infinite neutral grounded reservoir. Let us now write 
■0(r) = -00 + 5-0( r )- Notice that #0(r) is almost zero in the deep interior of the system 



Finite-Size Corrections for Coulomb Systems in the Debye-Hiickel Regime 



17 



and only has significant values near the boundaries. Let us suppose that the variations of 
5ip(r) are small, more precisely let us suppose that Si/j(r) = qg(nr) in two dimensions or 
5ip(r) = K,qg(K,r) in three dimensions with g some function of order one. Then expanding 
Si for small 5ip yields 

Si = lj 0V) + (fl&O'c) 0'(r) dr + 0(1*,) . (B.9) 

For the second part of the action the same expansion yields 

^ = " E ^ + ^ / A (^) + £~J ^(A - ^ + 0{T d d J (B.10) 

where we have know defined the inverse Debye length k in terms of the renormalized 
fugacities (* as k = a/ @ s d XL CaQa- Actually a closer inspection of equation (|B.10|) 
shows that the last term of S2 is actually of higher order that the two other terms. 
Indeed expanding Poisson-Boltzmann equation (jB.4|) for 8ip small shows that 



Therefore 



and 



JL 

2s d 



5^(A - k'W = 0(T% a ) (B.12) 



S * = -YlC a V + ^J v A(ty) + G(rJJ . (B.13) 

Notice now that only Si depends on <ft' and the result of the functional Gaussian 
integration over the field <$' will be the same as in section |21 equation (|15p. except 
that the fugacities £ a have to be replaced by £*. The term S2 does not depend on the 
field 0' and will give only an additive contribution to the grand potential. Finally we 
obtain for the grand potential 

n = kj f in (n - £) n e *) - ^ t e y c + ^ • (B.14) 

That is the same grand potential as before but evaluated for the fugacities C* instead 
of ( a plus a contribution 

S2 + J2CV = ~^ [ p{v)dv = -U Q Q. (B.15) 

We have used Poisson equation Atp = —s^p with p the average charge density of the 
system in the mean field approximation. The excess total charge is Q = f v p- Let us 
remark a few points about this term. The charge density p(r) is different from zero 
near the boundaries and a few Debye lengths away from the boudaries it vanishes. The 
charge Q is spread near the surface of the system. Therefore the additional contribution 
Qs to the grand potential is actually a surface contribution. This is even more clear if 
from equation (|B.13j) we use Gauss theorem to write Q,s a s 

^s = t^i V^(r)-rfS = ^/ a w {v)dS (B.16) 

^Sd JdV ^ JdV 



Finite-Size Corrections for Coulomb Systems in the Debye-Hiickel Regime 



18 



where a w = Sdd n ip is the surface charge density induced in the boundary walls of ideal 
conductor. This charge is external to the system. Since the ideal conductor is grounded 
and it is in total influence with the Coulomb system we have § dv a w dS = —Q, thus 
recovering equation (|B.15[) . This term could further be expressed as 

tts = ^- [ #(r)dr (B.17) 

2s d Jy 

where we have used equation (|B.11|) . Let us point out that this additional contribution 
Qs is n °t the naive electrostatic energy U e \ st = (1/2) f y p(r)?/>(r) dr. This can be 
checked for the two-dimensional two-component plasma near a planar wall made of 
ideal conductor. Using the results from Ref. [21], for small coulombic coupling, the 
mean field electric potential is ip(x) = ipo{l — e~ KX ), with x the distance from the wall. 
Then in this case = —Vts/2. 

In general in two or three dimensions Qs is of order Td, a - However for two 
dimensional systems when the pseudo-neutrality condition is satisfied, we have found a 

yf?rm 1 /2 

surface tension given by equation (jH7j) which is of lower order T 2 a - Therefore in two 
dimensions the surface contribution Q s to the grand potential due to the potential 
difference ipo between the system and the reservoir found in this appendix can be 
neglected in front of the surface tension given in equation (|H7|) . This fact has also been 
noticed in Ref. ^T7\ where the exact expression of the surface tension for a symmetric two- 
component plasma has been computed and its expansion for small coulombic coupling 
parameter shows that the potential ipo do not contribute in the dominant order. 

For two dimensional systems, we can conclude that the results of sections |2] 
and El also apply when the fugacities do not satisfy the pseudo-neutrality condition, 
provided that one replaces the original fugacities ( a with the renormalized fugacities 
C = C« e_/3<?a ^° which do satisfy the pseudo-neutrality condition. 

In three dimensions the situation is somehow different. By dimensional analysis 
one would expect that for a system of characteristic size R and with ip = the surface 
term of the grand potential will be proportional to k 2 R 2 , therefore of order r 3 Q ,.|| Then 
for a system with ip ^ the additional surface contribution Q$ found in this appendix 
will in principle contribute to the total surface tension (see, however, the footnote). 

Let us mention that if the potential difference is large one should go back to 
equations ()B.6j) and ()B.7j) and try to study the whole non-linear problem. There is an 
interesting regime where the fluctuations 0' around the mean field are small enough to 
expand the action S to the second order as S = Si + S2 but that the mean field ipo 
could be large and the further expansion of Si, equation (|B.9|) . and S 2 , equation (|B.13I) . 
is not possible. It is expected that some very interesting phenomena could occur in this 
non-linear regime, for instance renormalization and saturation of the surface charge Q 
and the potential i/j , like in the studies of highly charged colloids |IH1 H§1 loTIj. 

J Applying our method to three dimensional systems, actually gives a surface tension 7 = 
/3 _1 (k 2 /167t) ln(«;/A; max ) with /c max an ultraviolet cutoff. This result will be reported in a future work 
currently under preparation. Therefore it turns out that the dominant term in the surface tension is 
of order 1^3 Q lnr3 jQ and the correction Q$ is again of higher order. 
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To conclude this discussion we should point out a delicate point. The analysis 
done in this appendix is based on a mean field approximation: the function ip(r) and its 
constant value ipo inside the system are the solution of Poisson-Bolzmann equation (|B.4|I . 
In full generality they are different from the average electric potential inside the system. 
Only at the first order, in the ideal gas approximation (n a = £*) we can identify 
both. Nevertheless our goal was to obtain an estimation of the corrections to the 
grand-potential when the pseudo-neutrality condition is not satisfied, and based on 
this estimate we can conclude that these corrections are of higher order. 
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